#include <math.h>
#define NRANSI
#include "nrutil.h"
#include <stdio.h>

void fdjac(int n, double x[], double fvec[], double **df,
	void (*vecfunc)(int, double [], double []))
{
	int i,j;
	double h,temp,*f;
	extern int fdjac_counter;
    extern int verbose;
    extern int prodprob_in_jacobian;
    extern int usedcar_in_jacobian;
	extern double jaceps;

	fdjac_counter++;
	
	if(verbose) printf("\n\n\n*******\nstarting jacobian %d with jaceps %lf... ",fdjac_counter,jaceps);

	(*vecfunc)(n,x,fvec); //extra base eval in case passed fvec is not up to date
	
	//begin original fdjac code
	f=dvector(1,n);
	for (j=1;j<=n;j++) {
		temp=x[j];
		h=jaceps * fabs(temp);
		if (h == 0.0) h=jaceps;
		x[j]=temp+h;
		h=x[j]-temp;
		(*vecfunc)(n,x,f);
		x[j]=temp;
		for (i=1;i<=n;i++) df[i][j]=(f[i]-fvec[i])/h;
	}

	(*vecfunc)(n,x,fvec); //extra base eval to leave the most recent eval unchanged

	free_dvector(f,1,n);
	
}

#undef NRANSI
